############################################################################
########## Coup-Proofing: Latent Concept and Measurement
########## Replication, PSRM, Hwalmin Jin, 09/02/2023


######### 1) Computer specifications used to perform the Bayesian IRT estimation.
## Processor:	11th Gen Intel(R) Core(TM) i3-1115G4 @ 3.00GHz   3.00 GHz
## Installed RAM:	16.0 GB 
## System type:	64-bit operating system, x64-based processor
## Edition:	Windows 10 Enterprise
## Version:	21H2
## OS build:	19044.3086
## R version: 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Rstudio version: RStudio-2023.06.1-524


### Note (1): Regarding the replication

#Stan is designed to allow full reproducibility. However, this is only possible up to the external constraints imposed by floating point arithmetic.
#Stan results will only be exactly reproducible if all of the following components are identical:

#1) Stan version
#2) Stan interface (RStan, PyStan, CmdStan) and version, plus version of interface language (R, Python, shell)
#3) versions of included libraries (Boost and Eigen)
#4) operating system version
#5) computer hardware including CPU, motherboard and memory
#6) C++ compiler, including version, compiler flags, and linked libraries
#7) same configuration of call to Stan, including random seed, chain ID, initialization and data

#It doesn’t matter if you use a stable release version of Stan or the version with a particular Git hash tag. 
#The same goes for all of the interfaces, compilers, and so on. 
#The point is that if any of these moving parts changes in some way, floating point results may change.
#However, when reproduced with a system of a different specification, the replication results do not deviate by more than 5 percent from published results in the manuscript.

### Note (2): estimation time

#Based on the above computer specification, the estimation time is approximately 1 hour and 40 minutes per chain, 
#so the total time for four chains to complete the estimation is about 6 hours. 




######## 2) install rstan and preparation for IRT estimation.


install.packages("StanHeaders", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
install.packages("rstan", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))


## rstan version: 2.26.22
## StanHeaders version: 2.26.27

#example(stan_model, package = "rstan", run.dontrun = TRUE)


######### 3) Required R Packages

install.packages("rstantools")    ## version 2.3.1.1
install.packages("brms")          ## version 2.19.0
install.packages("reshape2")      ## version 1.4.4
install.packages("readr")         ## version 2.1.4
install.packages("dplyr")         ## version 1.1.2
install.packages("purrr")         ## version 1.0.1
install.packages("forcats")       ## version 1.0.0
install.packages("tidyr")         ## version 1.3.0
install.packages("modelr")        ## version 0.1.11
install.packages("ggdist")        ## version 3.3.0
install.packages("tidybayes")     ## version 3.0.4
install.packages("cowplot")       ## version 1.1.1
install.packages("ggrepel")       ## version 1.1.1
install.packages("gganimate")     ## version 1.0.8
install.packages("bayesplot")     ## version 1.10.0
install.packages("rstanarm")      ## version 2.21.4
install.packages("logr")          ## version 1.3.4
install.packages("rstudioapi")    ## version 0.15.0
install.packages("gridExtra")     ## version 2.3



library(rstan)
library(reshape2)
library(readr)
library(magrittr)
library(dplyr)
library(purrr)
library(forcats)
library(tidyr)
library(modelr)
library(ggdist)
library(tidybayes)
library(cowplot)
library(rstan)
library(ggrepel)
library(RColorBrewer)
library(gganimate)
library(bayesplot)
library(rstanarm)
library(ggplot2)
library(brms)
library(logr)
library(rstudioapi)
library(gridExtra)


### Generate log file using "logr"

## Create log file location
psrm_irt_logfile_R <- file.path("psrm_irt_logfile_R.log")

## Open log
lf <- log_open(psrm_irt_logfile_R)

log_open("psrm_irt_logfile_R.log")


## The below code will dump the entire program contents to the log. 
## The code lines will be prefixed with a right arrow (“>”) to distinguish these lines from the rest of your log.

log_code()

###You can create a log file for each command and data frame with the "put ( )" command.


# Get the directory of the current R script

put(dirname(rstudioapi::getSourceEditorContext()$path))

# Now you can perform your operations in the script's directory

put(getwd())

########### maximize print
put(options(max.print=3000))







######### 3. IRT Estimation


########## 1. IRT 2PL Model (Binary), based on "irt_coup_proofing_full.csv"
########## note: place a gamma prior (4, 3) for item discrimination parameter
########## saved file name: coup_proofing_rdata_psrm.RData



put(irt_coup_proofing_full.dat <- read.csv("irt_base_coup_full.csv",stringsAsFactors=F))


put(irt_coup_proofing_full_long <- melt(irt_coup_proofing_full.dat,id.vars=names(irt_coup_proofing_full.dat[,c(1:3)]),
                              measure.vars=names(irt_coup_proofing_full.dat[,4:ncol(irt_coup_proofing_full.dat)])))  ### convert to long format




######### IRT Model

put(formula_coup_proofing <- bf(
  value ~ beta + alpha * theta,
  theta ~ 0 + (1 | country_year),
  beta ~ 0 + variable,
  alpha ~ 0 + variable,
  nl = TRUE
))


put(prior_coup_proofing <-
  prior("normal(0, 10)", class = "b", nlpar = "beta", lb=-6, ub=6) +
  prior("gamma(4, 3)", class = "b", nlpar = "alpha", lb=0) +
  prior("normal(0, 1)", class = "sd", group = "country_year", nlpar = "theta"))




put(fit_coup_proofing <- brm(
  formula = formula_coup_proofing,
  data = irt_coup_proofing_full_long,
  family = brmsfamily("bernoulli", "logit"),
  prior = prior_coup_proofing,
  iter =5000,
))



######## Fitted Values 

put(fit_coup_proofing  %>% 
      spread_draws( r_country_year__theta[country_year,]) %>%
      head(10))

put(  fit_coup_proofing  %>% 
  spread_draws( r_country_year__theta[country_year,]) %>%
  head(10))

put(  fitted_coup_proofing <- spread_draws(fit_coup_proofing , r_country_year__theta[country_year,] ) %>%
  mutate(mean_coup_proofing = r_country_year__theta) %>%
  mean_qi(mean_coup_proofing))

put(write.csv(fitted_coup_proofing, "fitted_coup_proofing.csv" , row.names = FALSE))


#### Then, I imported the CSV file and converted it to the "fitted_coup_proofing.dta" format. 
#### And this file is merged into the stata replication file. 


########## Plotting Models   
put(posterior_fit_coup_proofing <- as.array(fit_coup_proofing))
put(dim(posterior_fit_coup_proofing))
put(dimnames(posterior_fit_coup_proofing))

put(color_scheme_set("blue"))

put(log_print(posterior_fit_coup_proofing))

########## 1) Item Easiness Parameter (Figure 1(a))
put(easiness_plot_coup_proofing <- mcmc_intervals(posterior_fit_coup_proofing, pars = c("b_beta_variablemil_mod_act_re", "b_beta_variablemil_int_act_re", 
                                                                              "b_beta_variablemil_moj_act_re", "b_beta_variablemil_mofa_act_re",
                                                                              "b_beta_variablemil_sec_other_act_re", "b_beta_variablemil_nonsec_act_re",
                                                                              "b_beta_variablepg", "b_beta_variableintel",
                                                                              "b_beta_variablepolice", "b_beta_variabletroops",
                                                                              "b_beta_variablemilitia", "b_beta_variablemilitary",
                                                                              "b_beta_variableborder", "b_beta_variablecw_mil",
                                                                              "b_beta_variablemilex_3yr_cut", "b_beta_variabletotal_purge_bi"), prob = FALSE,
                                             prob_outer = 0.95, # 95% intervals
                                             point_est = "mean", point_size = 2))



put(plot( easiness_plot_coup_proofing))



put(easiness_plot_coup_proofing + scale_y_discrete(labels = c("b_beta_variablemil_mod_act_re" ="Civilian Defense",
             "b_beta_variablemil_int_act_re" = "Civilian Interior",
             "b_beta_variablemil_moj_act_re" = "Civilian Justice",
             "b_beta_variablemil_mofa_act_re" = "Civilian Foreign Affairs",
             "b_beta_variablemil_sec_other_act_re" = "Civilian Security Position",
             "b_beta_variablemil_nonsec_act_re" = "Civilian Non-Security Position",
             "b_beta_variablepg"= "Presidential Guard",
             "b_beta_variableintel"= "Secret Police",
             "b_beta_variablepolice" = "Militarized Police",
             "b_beta_variabletroops" = "Interior Troops",
             "b_beta_variablemilitia" = "Militia",
             "b_beta_variablemilitary" = "Military Independent",
             "b_beta_variableborder"= "Borderguard",
             "b_beta_variablecw_mil"= "Military Counterweight",
             "b_beta_variablemilex_3yr_cut"= "Military Budget Reduction",
             "b_beta_variabletotal_purge_bi"= "Purge of Military Officers"),
  limits = c("b_beta_variablemil_mod_act_re",
             "b_beta_variablemil_int_act_re",
             "b_beta_variablemil_moj_act_re",
             "b_beta_variablemil_mofa_act_re",
             "b_beta_variablemil_sec_other_act_re",
             "b_beta_variablemil_nonsec_act_re",
             "b_beta_variablepg",
             "b_beta_variableintel",
             "b_beta_variablepolice",
             "b_beta_variabletroops",
             "b_beta_variablemilitia",
             "b_beta_variablemilitary",
             "b_beta_variableborder",
             "b_beta_variablecw_mil",
             "b_beta_variablemilex_3yr_cut",
             "b_beta_variabletotal_purge_bi")) +  scale_x_continuous("Estimate"))


####### To create a figure that looks like Figure 1 (a) in the manuscript, 

put(ggsave("easiness_static_figure1(a).pdf", height = 5.32, width = 7.8,  units = "in", scale = 0.7))


########## 2) Discrimination Parameter (Figure 1(b)) 


put(discrimination_plot_coup_proofing <- mcmc_intervals(posterior_fit_coup_proofing, pars = c("b_alpha_variablemil_mod_act_re", "b_alpha_variablemil_int_act_re", 
                                                                                    "b_alpha_variablemil_moj_act_re", "b_alpha_variablemil_mofa_act_re",
                                                                                    "b_alpha_variablemil_sec_other_act_re", "b_alpha_variablemil_nonsec_act_re",
                                                                                    "b_alpha_variablepg", "b_alpha_variableintel",
                                                                                    "b_alpha_variablepolice", "b_alpha_variabletroops",
                                                                                    "b_alpha_variablemilitia", "b_alpha_variablemilitary",
                                                                                    "b_alpha_variableborder", "b_alpha_variablecw_mil",
                                                                                    "b_alpha_variablemilex_3yr_cut", "b_alpha_variabletotal_purge_bi"), prob = FALSE,
                                                   prob_outer = 0.95, # 95% intervals
                                                   point_est = "mean",  point_size = 2) )




put(plot(discrimination_plot_coup_proofing))



put(discrimination_plot_coup_proofing + scale_y_discrete(
  labels = c("b_alpha_variablemil_mod_act_re" ="Civilian Defense",
             "b_alpha_variablemil_int_act_re" = "Civilian Interior",
             "b_alpha_variablemil_moj_act_re" = "Civilian Justice",
             "b_alpha_variablemil_mofa_act_re" = "Civilian Foreign Affairs",
             "b_alpha_variablemil_sec_other_act_re" = "Civilian Security Position",
             "b_alpha_variablemil_nonsec_act_re" = "Civilian Non-Security Position",
             "b_alpha_variablepg"= "Presidential Guard",
             "b_alpha_variableintel"= "Secret Police",
             "b_alpha_variablepolice" = "Militarized Police",
             "b_alpha_variabletroops" = "Interior Troops",
             "b_alpha_variablemilitia" = "Militia",
             "b_alpha_variablemilitary" = "Military Independent",
             "b_alpha_variableborder"= "Borderguard",
             "b_alpha_variablecw_mil"= "Military Counterweight",
             "b_alpha_variablemilex_3yr_cut"= "Military Budget Reduction",
             "b_alpha_variabletotal_purge_bi"= "Purge of Military Officers"),
  limits = c("b_alpha_variablemil_mod_act_re",
             "b_alpha_variablemil_int_act_re",
             "b_alpha_variablemil_moj_act_re",
             "b_alpha_variablemil_mofa_act_re",
             "b_alpha_variablemil_sec_other_act_re",
             "b_alpha_variablemil_nonsec_act_re",
             "b_alpha_variablepg",
             "b_alpha_variableintel",
             "b_alpha_variablepolice",
             "b_alpha_variabletroops",
             "b_alpha_variablemilitia",
             "b_alpha_variablemilitary",
             "b_alpha_variableborder",
             "b_alpha_variablecw_mil",
             "b_alpha_variablemilex_3yr_cut",
             "b_alpha_variabletotal_purge_bi")) +  scale_x_continuous("Estimate"))


####### To create a figure that looks like Figure 1 (b) in the manuscript, 

put(ggsave("discrimination_static_figure1(b).pdf", height = 5.32, width = 7.8,  units = "in", scale = 0.7))







######## Diagnostics ###########################################


put(prior_summary(fit_coup_proofing, all=FALSE))  


put(get_variables(fit_coup_proofing))


###### 1) Summary: Rhta, Bulk_ESS Tail_Ess
put(summary( fit_coup_proofing))





###### 2) distribution of rhats for the parameters (Figure 3 in Appendix)
put(rhats_fit_coup_proofing <- rhat(fit_coup_proofing))
put(print(rhats_fit_coup_proofing))

put(color_scheme_set("blue")) # see help("color_scheme_set")
put(mcmc_rhat(rhats_fit_coup_proofing))


####### To create a figure that looks like Figure 3 in Appendix, 

put(ggsave("Rhat_figure3appendix.pdf", height = 5.5, width = 7.8,  units = "in", scale = 0.7))


###### 3) MCMC Plots for the parameters (Figure 4 in Appendix)

put(mcmc_plot <- plot(fit_coup_proofing, ask=FALSE))

####### To create a figure that looks like Figure 4 in Appendix, 

put(Rplot_figure4appendix <- grid.arrange(grobs=mcmc_plot[2]))

put(ggsave("Rplot_figure4appendix.pdf", height = 5, width = 5.9,  units = "in", scale = 1.3, Rplot_figure4appendix))











########## 2. 2PL_2, Binary based on "irt_coup_proofing_full.csv"
########## note: place a normal prior (0, 10) for item discrimination parameter 
########## saved file name: coup_proofing_rdata_psrm_1.RData


put(irt_coup_proofing_full.dat <- read.csv("irt_base_coup_full.csv",stringsAsFactors=F))


put(irt_coup_proofing_full_long <- melt(irt_coup_proofing_full.dat,id.vars=names(irt_coup_proofing_full.dat[,c(1:3)]),
                                           measure.vars=names(irt_coup_proofing_full.dat[,4:ncol(irt_coup_proofing_full.dat)])))




######### IRT Model

put(formula_coup_proofing_1 <- bf(
  value ~ beta + alpha * theta,
  theta ~ 0 + (1 | country_year),
  beta ~ 0 + variable,
  alpha ~ 0 + variable,
  nl = TRUE
))


put(prior_coup_proofing_1 <-
  prior("normal(0, 10)", class = "b", nlpar = "beta", lb=-6, ub=6) +
  prior("normal(0, 10)", class = "b", nlpar = "alpha", lb=0) +
  prior("normal(0, 1)", class = "sd", group = "country_year", nlpar = "theta"))




put(fit_coup_proofing_1 <- brm(
  formula = formula_coup_proofing_1,
  data = irt_coup_proofing_full_long,
  family = brmsfamily("bernoulli", "logit"),
  prior = prior_coup_proofing_1,
  iter =5000,
))





put(prior_summary(fit_coup_proofing_1, all=FALSE))  

put(summary(fit_coup_proofing_1))


put(get_variables(fit_coup_proofing_1))



######## Fitted Values 

put(fit_coup_proofing_1  %>% 
  spread_draws( r_country_year__theta[country_year,]) %>%
  head(10))

put(fitted_coup_proofing_1 <- spread_draws(fit_coup_proofing_1 , r_country_year__theta[country_year,] ) %>%
  mutate(mean_coup_proofing_1 = r_country_year__theta) %>%
  mean_qi(mean_coup_proofing_1)) 


put(write.csv(fitted_coup_proofing, "fitted_coup_proofing_1.csv" , row.names = FALSE))

#### Then, I imported the CSV file and converted it to the "fitted_coup_proofing_1.dta" format.
#### And this file is merged into the stata replication file. 


########## Plotting Models   
put(posterior_fit_coup_proofing_1 <- as.array(fit_coup_proofing_1))
put(dim(posterior_fit_coup_proofing_1))
put(dimnames(posterior_fit_coup_proofing_1))


put(color_scheme_set("blue"))

########## 1) Item Easiness Parameter   
put(easiness_plot_coup_proofing_1 <- mcmc_intervals(posterior_fit_coup_proofing_1, pars = c("b_beta_variablemil_mod_act_re", "b_beta_variablemil_int_act_re", 
                                                                                                  "b_beta_variablemil_moj_act_re", "b_beta_variablemil_mofa_act_re",
                                                                                                  "b_beta_variablemil_sec_other_act_re", "b_beta_variablemil_nonsec_act_re",
                                                                                                  "b_beta_variablepg", "b_beta_variableintel",
                                                                                                  "b_beta_variablepolice", "b_beta_variabletroops",
                                                                                                  "b_beta_variablemilitia", "b_beta_variablemilitary",
                                                                                                  "b_beta_variableborder", "b_beta_variablecw_mil",
                                                                                                  "b_beta_variablemilex_3yr_cut", "b_beta_variabletotal_purge_bi"), prob = FALSE,
                                                     prob_outer = 0.95, # 95% intervals
                                                     point_est = "mean", point_size = 2))



put(plot( easiness_plot_coup_proofing_1))



put(easiness_plot_coup_proofing_1 + scale_y_discrete(labels = c("b_beta_variablemil_mod_act_re" ="Civilian Defense",
                                                                 "b_beta_variablemil_int_act_re" = "Civilian Interior",
                                                                 "b_beta_variablemil_moj_act_re" = "Civilian Justice",
                                                                 "b_beta_variablemil_mofa_act_re" = "Civilian Foreign Affairs",
                                                                 "b_beta_variablemil_sec_other_act_re" = "Civilian Security Position",
                                                                 "b_beta_variablemil_nonsec_act_re" = "Civilian Non-Security Position",
                                                                 "b_beta_variablepg"= "Presidential Guard",
                                                                 "b_beta_variableintel"= "Secret Police",
                                                                 "b_beta_variablepolice" = "Militarized Police",
                                                                 "b_beta_variabletroops" = "Interior Troops",
                                                                 "b_beta_variablemilitia" = "Militia",
                                                                 "b_beta_variablemilitary" = "Military Independent",
                                                                 "b_beta_variableborder"= "Borderguard",
                                                                 "b_beta_variablecw_mil"= "Military Counterweight",
                                                                 "b_beta_variablemilex_3yr_cut"= "Military Budget Reduction",
                                                                 "b_beta_variabletotal_purge_bi"= "Purge of Military Officers"),
                                                      limits = c("b_beta_variablemil_mod_act_re",
                                                                 "b_beta_variablemil_int_act_re",
                                                                 "b_beta_variablemil_moj_act_re",
                                                                 "b_beta_variablemil_mofa_act_re",
                                                                 "b_beta_variablemil_sec_other_act_re",
                                                                 "b_beta_variablemil_nonsec_act_re",
                                                                 "b_beta_variablepg",
                                                                 "b_beta_variableintel",
                                                                 "b_beta_variablepolice",
                                                                 "b_beta_variabletroops",
                                                                 "b_beta_variablemilitia",
                                                                 "b_beta_variablemilitary",
                                                                 "b_beta_variableborder",
                                                                 "b_beta_variablecw_mil",
                                                                 "b_beta_variablemilex_3yr_cut",
                                                                 "b_beta_variabletotal_purge_bi")) +  scale_x_continuous("Estimate"))




########## 2) Discrimination Parameter  


put(discrimination_plot_coup_proofing_1 <- mcmc_intervals(posterior_fit_coup_proofing_1, pars = c("b_alpha_variablemil_mod_act_re", "b_alpha_variablemil_int_act_re", 
                                                                                                        "b_alpha_variablemil_moj_act_re", "b_alpha_variablemil_mofa_act_re",
                                                                                                        "b_alpha_variablemil_sec_other_act_re", "b_alpha_variablemil_nonsec_act_re",
                                                                                                        "b_alpha_variablepg", "b_alpha_variableintel",
                                                                                                        "b_alpha_variablepolice", "b_alpha_variabletroops",
                                                                                                        "b_alpha_variablemilitia", "b_alpha_variablemilitary",
                                                                                                        "b_alpha_variableborder", "b_alpha_variablecw_mil",
                                                                                                        "b_alpha_variablemilex_3yr_cut", "b_alpha_variabletotal_purge_bi"), prob = FALSE,
                                                           prob_outer = 0.95, # 95% intervals
                                                           point_est = "mean",  point_size = 2) )




put(plot(discrimination_plot_coup_proofing_1))



put(discrimination_plot_coup_proofing_1 + scale_y_discrete(
  labels = c("b_alpha_variablemil_mod_act_re" ="Civilian Defense",
             "b_alpha_variablemil_int_act_re" = "Civilian Interior",
             "b_alpha_variablemil_moj_act_re" = "Civilian Justice",
             "b_alpha_variablemil_mofa_act_re" = "Civilian Foreign Affairs",
             "b_alpha_variablemil_sec_other_act_re" = "Civilian Security Position",
             "b_alpha_variablemil_nonsec_act_re" = "Civilian Non-Security Position",
             "b_alpha_variablepg"= "Presidential Guard",
             "b_alpha_variableintel"= "Secret Police",
             "b_alpha_variablepolice" = "Militarized Police",
             "b_alpha_variabletroops" = "Interior Troops",
             "b_alpha_variablemilitia" = "Militia",
             "b_alpha_variablemilitary" = "Military Independent",
             "b_alpha_variableborder"= "Borderguard",
             "b_alpha_variablecw_mil"= "Military Counterweight",
             "b_alpha_variablemilex_3yr_cut"= "Military Budget Reduction",
             "b_alpha_variabletotal_purge_bi"= "Purge of Military Officers"),
  limits = c("b_alpha_variablemil_mod_act_re",
             "b_alpha_variablemil_int_act_re",
             "b_alpha_variablemil_moj_act_re",
             "b_alpha_variablemil_mofa_act_re",
             "b_alpha_variablemil_sec_other_act_re",
             "b_alpha_variablemil_nonsec_act_re",
             "b_alpha_variablepg",
             "b_alpha_variableintel",
             "b_alpha_variablepolice",
             "b_alpha_variabletroops",
             "b_alpha_variablemilitia",
             "b_alpha_variablemilitary",
             "b_alpha_variableborder",
             "b_alpha_variablecw_mil",
             "b_alpha_variablemilex_3yr_cut",
             "b_alpha_variabletotal_purge_bi")) +  scale_x_continuous("Estimate"))



####### export size: pdf (5.32 x 3.85)



######## Diagnostics 



###### 1) Summary: Rhta, Bulk_ESS Tail_Ess
put(summary( fit_coup_proofing_1))




####### 2) MCMC Plots for the parameters
#mcmc_trace(posterior_fit_coup_proofing_effort, pars = c("b_eta_variableMilitary_Defense", "b_eta_variableMilitary_Interior"))

put(plot(fit_coup_proofing_1, ask=FALSE))



######## 3) distribution of rhats for the parameters
put(rhats_fit_coup_proofing_1 <- rhat(fit_coup_proofing_1))
put(print(rhats_fit_coup_proofing_1))

put(color_scheme_set("blue")) # see help("color_scheme_set")
put(mcmc_rhat(rhats_fit_coup_proofing_1))







########## 3. 2PL_2, Binary based on "irt_coup_proofing_effort_full_exclude.csv"
########## Note: exclude civilian elites in foreign affairs position, justice, and non-security cabient posts
########## saved file name: coup_proofing_rdata_psrm_exclude.RData


put(irt_coup_proofing_full_exclude.dat <- read.csv("irt_base_coup_full_exclude.csv",stringsAsFactors=F))


put(irt_coup_proofing_full_exclude_long <- melt(irt_coup_proofing_full_exclude.dat,id.vars=names(irt_coup_proofing_full_exclude.dat[,c(1:3)]),
                                           measure.vars=names(irt_coup_proofing_full_exclude.dat[,4:ncol(irt_coup_proofing_full_exclude.dat)])))




######### IRT Model

put(formula_coup_proofing_exclude <- bf(
  value ~ beta + alpha * theta,
  theta ~ 0 + (1 | country_year),
  beta ~ 0 + variable,
  alpha ~ 0 + variable,
  nl = TRUE
))


put(prior_coup_proofing_exclude <-
  prior("normal(0, 10)", class = "b", nlpar = "beta", lb=-6, ub=6) +
  prior("gamma(4, 3)", class = "b", nlpar = "alpha", lb=0) +
  prior("normal(0, 1)", class = "sd", group = "country_year", nlpar = "theta"))




put(fit_coup_proofing_exclude <- brm(
  formula = formula_coup_proofing_exclude,
  data = irt_coup_proofing_full_exclude_long,
  family = brmsfamily("bernoulli", "logit"),
  prior = prior_coup_proofing_exclude,
  iter =5000,
))





put(prior_summary(fit_coup_proofing_exclude, all=FALSE))  

put(summary(fit_coup_proofing_exclude))


put(get_variables(fit_coup_proofing_exclude))



######## Fitted Values 

put(fit_coup_proofing_exclude  %>% 
  spread_draws( r_country_year__theta[country_year,]) %>%
  head(10))

put(fitted_coup_proofing_exclude <- spread_draws(fit_coup_proofing_exclude , r_country_year__theta[country_year,] ) %>%
  mutate(mean_coup_proofing_exclude = r_country_year__theta) %>%
  mean_qi(mean_coup_proofing_exclude)) 


put(write.csv(fitted_coup_proofing, "fitted_coup_proofing_exclude.csv" , row.names = FALSE))

#### Then, I imported the CSV file and converted it to the "fitted_coup_proofing_exclude.dta" format.
#### And this file is merged into the stata replication file. 



########## Plotting Models   
put(posterior_fit_coup_proofing_exclude <- as.array(fit_coup_proofing_exclude))
put(dim(posterior_fit_coup_proofing_exclude))
put(dimnames(posterior_fit_coup_proofing_exclude))


put(color_scheme_set("blue"))

########## 1) Item Easiness Parameter (Figure 1(a) in Appendix)   
put(easiness_plot_coup_proofing_exclude <- mcmc_intervals(posterior_fit_coup_proofing_exclude, pars = c("b_beta_variablemil_mod_act_re", "b_beta_variablemil_int_act_re", 
                                                                                                  "b_beta_variablemil_sec_other_act_re", 
                                                                                                  "b_beta_variablepg", "b_beta_variableintel",
                                                                                                  "b_beta_variablepolice", "b_beta_variabletroops",
                                                                                                  "b_beta_variablemilitia", "b_beta_variablemilitary",
                                                                                                  "b_beta_variableborder", "b_beta_variablecw_mil",
                                                                                                  "b_beta_variablemilex_3yr_cut", "b_beta_variabletotal_purge_bi"), prob = FALSE,
                                                     prob_outer = 0.95, # 95% intervals
                                                     point_est = "mean", point_size = 2))



put(plot( easiness_plot_coup_proofing_exclude))



put(easiness_plot_coup_proofing_exclude + scale_y_discrete(labels = c("b_beta_variablemil_mod_act_re" ="Civilian Defense",
                                                                 "b_beta_variablemil_int_act_re" = "Civilian Interior",
                                                                 "b_beta_variablemil_sec_other_act_re" = "Civilian Security Position",
                                                                 "b_beta_variablepg"= "Presidential Guard",
                                                                 "b_beta_variableintel"= "Secret Police",
                                                                 "b_beta_variablepolice" = "Militarized Police",
                                                                 "b_beta_variabletroops" = "Interior Troops",
                                                                 "b_beta_variablemilitia" = "Militia",
                                                                 "b_beta_variablemilitary" = "Military Independent",
                                                                 "b_beta_variableborder"= "Borderguard",
                                                                 "b_beta_variablecw_mil"= "Military Counterweight",
                                                                 "b_beta_variablemilex_3yr_cut"= "Military Budget Reduction",
                                                                 "b_beta_variabletotal_purge_bi"= "Purge of Military Officers"),
                                                      limits = c("b_beta_variablemil_mod_act_re",
                                                                 "b_beta_variablemil_int_act_re",
                                                                 "b_beta_variablemil_sec_other_act_re",
                                                                 "b_beta_variablepg",
                                                                 "b_beta_variableintel",
                                                                 "b_beta_variablepolice",
                                                                 "b_beta_variabletroops",
                                                                 "b_beta_variablemilitia",
                                                                 "b_beta_variablemilitary",
                                                                 "b_beta_variableborder",
                                                                 "b_beta_variablecw_mil",
                                                                 "b_beta_variablemilex_3yr_cut",
                                                                 "b_beta_variabletotal_purge_bi")) +  scale_x_continuous("Estimate"))


####### To create a figure that looks like Figure 1 (a) in Appendix, 

put(ggsave("easiness_exclude_static_figure1(a)appendix.pdf", height = 5.32, width = 7.8,  units = "in", scale = 0.7))


########## 2) Discrimination Parameter(Figure 1 (b) in Appendix)


put(discrimination_plot_coup_proofing_exclude <- mcmc_intervals(posterior_fit_coup_proofing_exclude, pars = c("b_alpha_variablemil_mod_act_re", "b_alpha_variablemil_int_act_re", 
                                                                                                        "b_alpha_variablemil_sec_other_act_re", 
                                                                                                        "b_alpha_variablepg", "b_alpha_variableintel",
                                                                                                        "b_alpha_variablepolice", "b_alpha_variabletroops",
                                                                                                        "b_alpha_variablemilitia", "b_alpha_variablemilitary",
                                                                                                        "b_alpha_variableborder", "b_alpha_variablecw_mil",
                                                                                                        "b_alpha_variablemilex_3yr_cut", "b_alpha_variabletotal_purge_bi"), prob = FALSE,
                                                           prob_outer = 0.95, # 95% intervals
                                                           point_est = "mean",  point_size = 2) )




put(plot(discrimination_plot_coup_proofing_exclude))



put(discrimination_plot_coup_proofing_exclude + scale_y_discrete(
  labels = c("b_alpha_variablemil_mod_act_re" ="Civilian Defense",
             "b_alpha_variablemil_int_act_re" = "Civilian Interior",
             "b_alpha_variablemil_sec_other_act_re" = "Civilian Security Position",
             "b_alpha_variablepg"= "Presidential Guard",
             "b_alpha_variableintel"= "Secret Police",
             "b_alpha_variablepolice" = "Militarized Police",
             "b_alpha_variabletroops" = "Interior Troops",
             "b_alpha_variablemilitia" = "Militia",
             "b_alpha_variablemilitary" = "Military Independent",
             "b_alpha_variableborder"= "Borderguard",
             "b_alpha_variablecw_mil"= "Military Counterweight",
             "b_alpha_variablemilex_3yr_cut"= "Military Budget Reduction",
             "b_alpha_variabletotal_purge_bi"= "Purge of Military Officers"),
  limits = c("b_alpha_variablemil_mod_act_re",
             "b_alpha_variablemil_int_act_re",
             "b_alpha_variablemil_sec_other_act_re",
             "b_alpha_variablepg",
             "b_alpha_variableintel",
             "b_alpha_variablepolice",
             "b_alpha_variabletroops",
             "b_alpha_variablemilitia",
             "b_alpha_variablemilitary",
             "b_alpha_variableborder",
             "b_alpha_variablecw_mil",
             "b_alpha_variablemilex_3yr_cut",
             "b_alpha_variabletotal_purge_bi")) +  scale_x_continuous("Estimate"))


####### To create a figure that looks like Figure 1 (b) in Appendix, 

put(ggsave("discrimination_exclude_static_figure1(b)appendix.pdf", height = 5.32, width = 7.8,  units = "in", scale = 0.7))


######## Diagnostics 



###### 1) Summary: Rhta, Bulk_ESS Tail_Ess
put(summary(fit_coup_proofing_exclude))




####### 2) MCMC Plots for the parameters
#mcmc_trace(posterior_fit_coup_proofing_effort, pars = c("b_eta_variableMilitary_Defense", "b_eta_variableMilitary_Interior"))





######## 3) distribution of rhats for the parameters
put(rhats_fit_coup_proofing_exclude <- rhat(fit_coup_proofing_exclude))
put(print(rhats_fit_coup_proofing_exclude))

put(color_scheme_set("blue")) # see help("color_scheme_set")
put(mcmc_rhat(rhats_fit_coup_proofing_exclude))






## Close log
log_close()

## View results
writeLines(readLines(lf))


   